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SUMMARY 


The flow over a helicopter rotor blade in forward flight is an important 
example of three-dimensional tine-dependent flow. The boundary layers on the 
rotor blade set loss levels and control retreating blade stall. As a conse- 
quence there is considerable interest in developing a numerical scheme for 
solving the time-dependent viscous compressible three-dimensional flow to aid 
in the design of helicopter rotors. In the present report candidate numerical 
algorithms are examined to determine their overall suitability for the 
efficient and routine solution of an appropriate system of partial differential 
equations. It is concluded that a consistently split time-linearized block 
implicit scheme using either quintic B-spline collocation or the generalized 
operator compact implicit approach to generate a fourth order accurate 
algorithm is particularly well suited for use on the present problem. High 
cell Reynolds number behavior leads to favoring the generalized operator 
compact implicit approach over the quintic B-spline collocation method. 



INTRODUCTION 



The behavior of boundary layers on wings and bodies has long been of 
interest to aerodynamicists. In both steady and unsteady flows the boundary 
layers are known to govern a major portion of the losses and to significantly 
Influence the vehicle lift and moment coefficients. When the flow is steady, 
boundary layer prediction schemes based on numerical solution to the governing 
partial differential equations of motion have reached a high level of sophis- 
tication and predictive accuracy, even in three space dimensions. In 
unsteady flows, ouch as are commonly encountered in rotary winged aircraft, 
some progress has been made in two space dimensions but little to date has 
appeared on unsteady three-dimensional boundary layers. 

Two particular problems arise with time-dependent three-dimensional 
boundary layers relative to the steady case. The first of these is the 
rather obvious one of time integration with its added requirements of transient 
accuracy coupled with an increase in the computational labor. The second of 
these is the so-called negative cross flow problem, which to some extent has 
troubled the steady boundary layer prediction schemes. Kendall et al., 

(Ref. 1) discuss the negative cross flow problem for steady three-dimensional 
boundary layers in a very illuminating fashion. This particular problem 
arises when the spanwise component of velocity changes sign and will be 
discussed in detail subsequently. Because of the interest by external 
aerodynamlclsts in swept wing boundary layers where the negative cross flow 
problem (in this case flow from tip to root) is not usually encountered, the 
negative cross flow problem has not received a great deal of attention to 
date. However in transient flows, particularly those encountered on rotor 
blades in forward flight, negative cross flows are frequently encountered. 

For instance, the advancing rotor blade has cross flows of one sign during 
the first ninety degrees of rotation and these can change sign over part of 
the blade during the second ninety degrees. 

Thus to be of practical value, time-dependent three-dimensional boundary 
layer prediction schemes require high computational efficiency and transient 
accuracy coupled to the ability to treat arbitrary cross flow profiles. 

These attributes arc not available in any existing available computer code 
and hence in view of the potential use for a code of this type its development 
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is appropriate and timely. In this report the initial phase in the development 
of an efficient time-dependent three-dimensional boundary layer prediction 
procedure is investigated, namely, the choice of the computational algorithm 
and spatial differencing technique. 
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Background 

In this section we discuss the requirements of a three-dimensional 
unsteady boundary layer code and demonstrate how the physics of the flow 
influences the choice of a computational technique. 

Three-dimensional boundary layers occur on the wings and fuselages of 
both conventional and rotary wing aircraft. In both types of vehicles, the 
boundary layers are important in setting loss levels and determining useful 
operating ranges. As is well known, boundary layers are sensitive to pressure 
gradients. In time-dependent flow the temporal acceleration terms appear in 
the momentum equation in a form very similar to the conventional imposed 
pressure gradient and so for qualitative evaluation purposes can be regarded 
as ’pseudo' or ’auxiliary’ pressure gradients. Viewed in this manner the 
temporal acceleration terms are likely to be able to influence quantities of 
practical importance such as skin friction, displacement thickness and the 
onset of separation. At the range of frequencies typically encountered in 
rotary wing aircraft aerodynamic problems, it is clear, for instance, from 
the very thorough review of McCroskey (Ref. 2), that very significant transient 
boundary layer effects can be observed. 

In examining the flow problems of practical interest such as loss levels 
or the onset of separation it is evident that all three space dimensions must 
be considered. In conventional aircraft the sweep effect is of interest and 
inherently three-dimensional. In rotary wing aircraft in forward flight 
clearly very substantial transient changes occur in what might be termed the 
local sweep angle. However generally speaking, the boundary layers remain 
thin unless catastrophic flow separation occurs or the flow at the wing or 
rotor tip is considered. As a consequence it might be supposed that the usual 
three-dimensional thin boundary sheet approximations (Nash and Patel, Ref. 3) 
could be used to produce a valid set of governing equations. Fortunately some 
improvements in thin boundary sheet approximations are possible ae a result 
of having to treat the negative cross flow problem mentioned earlier. 

The negative cross flow problem is best explained in a somewhat intuitive 
manner, and for steady boundary layers a very good physical description of the 
problem is given by Kendall et al., (Ref. 1). Looking at the suction surface of 
a conventional swept back wing the boundary layer cross flow, w, is usually 
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outward in the z positive direction along the span f’-om root to tip. Thus 
conventional steady boundary layer integration schemes have developed by 
forward marching the streamwise velocity u in the streamwise x direction and 
simultaneously marching out along the span in the z positive direction. In 
view of the physics of the problem, the spanwise marching scheme does not 
normally encounter negative w, i.e., spanwise inflow. This i3 very fortunate 
because it is difficult, indeed it could be argued impossible, to structure a 
physically satisfactory unconditionally stable scheme which permits forward 
marching in the spanwise direction with a negative w cross flow. At least 
Intuitively the problem of negative cross flow implies information being 
transferred upstream against the spanwise marching direction. Conventional 
stability analyses confirm the inability to forward march into regions of 
significant negative w. From experience with attempts to march the two- 
dimensional boundary layer equations into a region of separated flow and its 
obvious relationship to the negative cross flow problem, it is not surprising 
that spanwise marching into a negative cross flow region is not accomplished 
without special treatment, for instance the Krause "zig-zag" scheme (Ref. 4). 
Recently conventional boundary layer developers have been turning to a span- 
wise as well as normal implicit formulation to remove the restriction of only 
positive cross flows (Kendall et al., Ref. 1). With a spanwise implicit 
formulation spanwise diffusion is allowed, and the resulting implicit system 
of equations can be treated by direct elimination (Ref. 1), by a predictor- 
corrector iterative approach (Ref. 5), or by the process of matrix splitting which 
reduces the matrix elimination labor (Refs. 6 and 7). Lin and Rubin (Ref. 5) 
in their predictor-corrector boundary region solutions for flow over a yawed 
cone at moderate incidence showed that allowing diffusion in the spanwise 
direction not only eliminates the problems associated with negative cross flow, 
but improves upon the solutions obtained by three-dimensional boundary layer 
techniques. Again intuitively a spanwise implicit construction permits 
information transfer in either direction. Boundary conditions applied at the 
tip can influence the flow inboard, if required by the physics of the flow. 

For these reasons the implicit spanwise construction has been a feature of 
the three-dimensional duct flow analysis of Briley (Ref. 6} and McDonald and 
Briley (Ref. 7). Based on the experience in Refs. 6 and 7, the additional 
computational effort resulting from a spanwise implicit formulation could be 
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as little as a 20Z increase relative to the explicit spanwise inarching approach. 
The extension of the steady three-dimensional boundary layer equations to allow 
spanwise diffusion is easily accomplished, and in view of the improved physical 
representation which thus follows, it is recommended. As a matter of course it 
has been assumed that normal to the wall an implicit formulation would be 
structured. In recent years for boundary layer type problems there has been 
little dispute as to the efficiency gains to be had from an implicit formula- 
tion normal to the wall (Ref. 8). However in the streamwise direction for 
steady two-dimensional flow, the equations are normally forward marched and 
the implicit stability obtained entirely from being implicit in the normal to 
the wall direction. 

For unsteady boundary layers the problem can again be formulated in 
either an explicit or implicit manner. As with spatial marching of steady 
boundary layers for usual aerodynamic applications, the locally refined 
spatial mesh required to define a (turbulent) boundary layer influenced by a 
(transiently) varying pressure distribution, when an explicit (stability 
restricted) scheme is employed, results in a maximum time step that is much 
less than the time scale of the physical processes of interest. Thus for 
solving unsteady boundary layers of the type usually encountered in rotary 
winged aircraft an implicit formulation is desirable. Since in time-dependent 
flow diffusion in the streamwise direction is normally negligible due to the 
usual boundary layer approximations, it is possible to formulate an implicit 
time-dependent scheme that retains the implicit structure in the spanwise 
and normal directions (which was found desirable for the steady boundary 
layer) and march the solution in the streamwise direction. 

As mentioned earlier the streamwise marching sweep would probably require 
less computational effort by about 20% than a fully implicit formulation and 
of course less storage. However since the solution is being time marched, the 
opportunity to use a streamwise implicit formulation at roughly the same cost 
as the streamwise marching sweep does arise. If one does perform a streamwise 
marching sweep, then the linearization of nonlinear terms is performed about 
the known spatial marching level. If a fully implicit structure is adopted, 
then full time linearization can be utilized. That is the linearization of the 
nonlinear terms is performed about the known time level. As is pointed out in 
Ref. 7, it is easier to obtain a consistent spatial-temporal order accurate 
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linearization by marching in time than in space (in time the nonlinear marching 
derivatives have the form Pu^ whereas in space marching they have the form 
pUjUj). Further by structuring implicitly in the space marching direction, 
(small) regions of axial reverse flow would be permitted. As a result of 
these combined benefits of linearization and separation, a fully implicit 
structure is advocated. 

Transient calculations mean that, in essence, a full three-dimensional 
spatial integration is carried out at each time step. Thus, spatial accuracy 
is very important to minimize the spatial grid point density for efficiency 
since many time steps are contemplated in a given cycle. In order to get the 
most out of a given spatial difference formula, the errors from representing 
nonlinear terms by linear combinations of terms should be less than or equal to 
the spatial discretization errors. If the linearization introduces a greater 
error than the spatial differencing, then either a coarser spatial mesh could 
be used, or iteration, or some form of linearization improvement is called 
for. Iteration across a time step is not recommended since this only reduces 
the linearization error and computationally costs as much as a complete time 
step. Cutting back the time step would be preferable to iterating to preserve 
the linearization error at some acceptable level, since cutting back on the 
time step would improve both ' the transient error and the linearization error. 

To obtain a linearization, which introduces errors of at most the same as the 
spatial difference formulae, a Taylor series expansion about the known time 
level can be performed. This process clearly demands a formal block, i.e. , 
coupled, treatment of the system of equations. For instance in the streamwise 
momentum equation a typical term is linearized: 

, .n+i n+i n . n n+i n n .,2, 

( uw) = u w + u w -uw+OlAt) 

and clearly one cannot lag w n+ ^ at the old time level n without introducing a 
first order time error in order to get an uncoupled system, i.e., w 11 " 1 ^ not 
appearing in the streamwise momentum equation. Thus formal linearization and 
consideration of the resulting errors indicate the coupled system ought to be 
treated from the accuracy point of view. This is further reinforced when it 
is realized that direct elimination of block, i.e., coupled, banded systems 
are not computationally expensive compared to the iterative solution of an 
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uncoupled system. The linearization technique is described In detail in 
Appendix A. 

Additionally a second type of approximation arises unconnected with 
linearization but arising from basic coupling terms in the original equations 
and if indeed some terms in an equation are time lagged in order to uncouple 
the equation system and these terms are of equal importance to the terms 
retained, then again an iterative updating is called for in order to achieve 
stability, accuracy and consistency. (This could be termed ad hoc equation 
uncoupling). Blottncr (Ref. 8) has shown that many iterations around the 
ad hoc uncoupled set (>10) is sometimes required in order to achieve an overall 
solution accuracy commensurate with the local difference molecule accuracy. 
Thompson and MacDonald (Ref. 9) found in a three-dimensional momentum integral 
procedure that a lagged sequential iterative type of calculation would not, 
in a number of instances, even converge. The weight of opinion definitely 
favors the block coupled approach. 

As a general observation, care is required to obtain acceptable transient 
accuracy for long time integration with conventional finite difference schemes. 
A Crank-I.'icolson centered time implicit scheme for instance, although second 
order in time, shows quite a dispersion problem (relative to other schemes) on 
the simple pure convection problem. However the problem of transient accuracy 
is significantly reduced in the typical boundary layer problem since the time 
dependency is continuously input through initial and boundary conditions and 
relatively the concern is with 'short* time integrations. The computational 
problem is more of uhat the phase lag of the wall shear is, relative to the 
prescribed free stream disturbance, than concern over the convection velocity 
of a wave in a shear after a long propagation tine. The interest is in forced 
oscillations with a minimum scale of the boundary layer thickness over a few 
cycles of the motion, just enough to obtain repetition cyclically. It is 
therefore expected that a significant dispersion problem will not arise with 
a conventional implicit scheme. 

The equation system which will be considered is formally of block size 
four, consisting of the continuity and two momentum equations and an 'energy' 
equation for p, u, v, w (p is specified everywhere). If constant stagnation 
temperature is assumed, p, u and w are related by an algebraic equation and 
the problem C3n be reduced to a block-three system rather than block-four upon 
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option, with significant reduction in computer tine. Matrix splitting 
techniques such as those described in Refs. 6 and 7 have considerable potential 
to reduce the computational labor of solving the block implicit system of 
(linear) algebraic equations which result from discretizing the governing 
equations. Schemes of this general type are termed split linearized block 
implicit or for brevity split LB1 schemes and are reviewed in detail by 
Briley and McDonald (Ref. 10). Ulth a careful ordering of the sweeps with 
a split LBI scheme it is possible to use a block-two with the third equation 
uncoupled on two of the sweeps so that this would be a potential major 
advantage of a split LBI approach. 

The ultimate goal of the "optimum" scheme is to diminish both storage 
requirements and running times in order to achieve a desired accuracy level. 
Although schemes can be constructed to satisfy cither one or both of these 
goals, the robustness of a method can only be verified by considering its 
applicability to a general class of problems. Our concern here is with an 
approximate form of the three-dimensional unsteady compressible Navier-Stokes 
equations, so that the method chosen will by necessity be required to treat a 
coupled system of nonlinear partial differential equations. At the outset the 
following observations can be made concerning the characteristics of the method. 

1. Implicit methods (preferably in all three spatial directions) are desired 
in order to eliminate stability restrictions and permit solutions with both 
positive and negative streamwir.e and spanwise velocities. 

2. The nature of the nonlinear coupling of the variables in the governing 
equations require that the equations be solved coupled. 

3. Iteration should be avoided and "time linearization" procedures employed 
(a discussion of this point is given in the following section). 

A. The method should allow for general boundary conditions. 

5. The method should allow the flexibility of incorporating higher order 

spatial differencing methods. 

Although Item 5 pertains to spatial differencing, in the next section we 
will demonstrate that it also plays a crucial role in choosing the type of 
temporal scheme. It would thus appear that first an overall temporal discreti- 
zation procedure must be chosen and only then an efficient spatial scheme that 
could be incorporated with it. 

The arguments above, linearization, stability considerations and physics 
lead us to suggest that the current state of the art dictates that the recommended 
scheme should be in the framework of a tine linearized block implicit method. 
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With the foregoing as background indicating the recommended overall 
approach to the problem of developing an efficient accurate computer code to 
predict time-dependent viscous flows, attention is turned to the specific prob- 
lem dealt with in this report. This report concerns itself with first which 
overall numerical algorithm to adopt and secondly which spatial differencing 
scheme to employ or what basis of approximating solution functions to use. 

In the next section we will discuss three methods for solving general 
systems of multidimensional parabolic equations, the consistently split block 
implicit scheme, Rubin's Predictor-Corrector technique and the Hopscotch 
algorithm. The following two sections will describe two spatial differencing 
methods, which we classify as Q-R operator schemes and basis function schemes 
(which deal mainly with B-splineo). Applications of these methods to model 
two point boundary value problems and to a coupled system of two nonlinear one- 
dimensional parabolic equations possessing a "three-dimensional boundary 
layer-like behavior" are then given. 

Temporal Schemes 

Consider a system of three-dimensional nonlinear parabolic differential 
equations 

<£, m %(<}>, <P x ,4\ x ,4>y,<i>y y ,<t> z ,4> zz ,*, y,z,t) (1) 

where ^ is a vector of unknowns ($j, . With the equations appropriately 

linearized, at each tine step, a system of N - 3(1-1) (J-l) (K-l) linear equations 
result, where I, J, K are the number of intervals In the x, y, z directions, 
respectively. Direct inversion of the system is not practical in three space 
dimensions except for extremely coarse meshes since the operation count is 
proportional to N /3 . Although higher order spatial schemes allows one to 

reduce the number of grid points, results of model problems indicate that for 
the range of accuracy desired, one can expect at best a reduction of a factor 
of four in grid points in each coordinate direction. As significant as this may 
seem, it nay not appreciably affect the overall computation unless the original 
problem is reduced to a more tractable form. 

*For I**J”K"10 N° 2187 and the matrix inversion operation count is 

proportional to 10^®. 

For I ■* J - K • 50 N *= 352947 and the matrix inversion operation count is 

proportional to 10* 6 . 
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The usual procedure lo to transform the original problem into a sequence 
of simpler problems or into a sequence of one-dimensional problems. We 
consider the linearized version of Eq. (1) (see Appendix A for a description of 
the time linearization procedure) with the differential operators identified 
with their coordinate direction, namely 

{^ n+, -^ n )/At^; + ^ n y+ ^;)(/?/ 4i + (i-/3)^ n ) +s ft ( 2 ) 

where superscripts indicate the time level, i.e., t n ■ nAt, S n is a source 
term and 0 5 0 £ 1 is a parameter allowing one to center the time step, i.e., 

6 ■ 0 corresponds to a forward difference, 0 “ 1/2 to Crank-HicolGon and 
0 ■ 1 to a backward difference. 

Consistently Split Block Implicit Scheme 

We will first consider the consistently split block Implicit scheme. 

Solution of Eq. (2) is accomplished by application of a generalization to 
systems of PDE's of an alternating-direction implicit (ADI) technique for 
parabolic-hyperbolic equations. The original ADI method was introduced by 
Peaceman and Rachford (Ref. 11) and Douglas (Ref. 12); however, the 
alternating-direction concept has since been expanded and generalized. A 
discussion of various alternating-direction techniques is given by Mitchell 
(Ref. 13), Yanenko (Ref. 14) and more recently by Briley and McDonald (Ref. 10). 

The present technique is simply an application of a generalization of 
the procedure developed by Douglas and Gunn (Ref. 15) for generating consis- 
tently split ADI schemes as perturbations of fundamental implicit difference 
schemes such as the backward-difference or Crank-Nicolson schemes in its 
natural extension to systems of partial differential equations. In this 
context a consistent scheme is one where the intermediate levels represent 
a discrete approximation to the governing equations whose truncation error can 
be made to vanish as the time or spatial mesh is arbitrarily reduced. Consistency 
in this sense is a very valuable property as it can greatly simplify the 
accurate implementation of boundary conditions (Ref. 10). 

For the present, it will be assumed that _^($) contains derivatives of first 
and second order with respect to the coordinate direction, but no mixed deriva- 
tives, Mixed derivatives are allowable within the formal framework but unless 
they are important they arc best treated explicitly (lagging) or by extrapolation. 
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The Douglas-Gunn representation of Eq. (2) can be written as the following 
three-step solution procedure: 

(<£*-<£ n )/A t + [l\-(3)\ + 2> y +% z ]<£ n + s" 

( «/»"- <£ n ) /At • #& x f +/32> y <£ # * + [ ( I - pn\ + 2> y )+2> 2 ] 4? +s n 

( <£*" - 4?) / At • +p2> y <T + + [ ( i - £>a x + ^ y + 2> z )] /+s n 

<£ n+l . ^•••+o(At 3 ) 

(3) 


In increment form Eq. (3) reduces to the algorithm given by Briley and 
McDonald (Refs. 10 and 17) for solving the compressible time-dependent three- 
dimensional Navier-Stokes equations, viz., 

( I -/3At2> x )(<P* - <£ n ) - At(S> x +2> y +\)4> n + AtS n 


% 


i 


t * 


(I -/3AO> y )(<£** -<£")- <£* - <p n 


(4) 


«£"*'.<£***+ 0 (At 3 ) 

A ** 

where $ and <J> are intermediate solutions. Each of Eqs. (4) can be written 

in narrow block-banded matrix form and solved by efficient block-elimination 
* ** 

methods. If $ and <{> are eliminated, Eqs. (4) become 

(i - £At2> x )(i - /3At2> y )( i - /?At2> z )(<£ n+ ' - *"> - At [c& x + 2> y +\)4> n + s n ] 

(5) 
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If the multiplication on the left-hand side of Eq. (5) is performed, it 
becomes apparent that Eq. (5) approximates Eq. (2) to order (At) . Although 
the stability of Eqo. (4) has not been established in circumstances suffi- 
ciently general to encompass the Navier-Stokes equations, it is often suggested 
(e.g. , Richtmyer and Morton, Ref. 16, p. 215) that the scheme is stable and 
accurate under conditions more general than those for which rigorous proofs 
are available. This latter notion was adopted here as a working hypothesis 
supported by favorable results obtained in actual computations (e.g., Refs. 

17, 18 and 19). 

Several observations can be made concerning Eqs. (3) and (4). 

1. The system of three-dimensional equations has been reduced to three 
systems of one-dimensional equations. 

2. The inversion of the total system is now at most approximately 108 

(IJ + JK + IK) operations compared to approximately (3IJK) J operations for the 
direct Inversion problem. 

3. The first step involves at least 402 - 50% of the operations (Ref. 7). 

4. The method does not have a CFL stability condition and for B 1 1/2 is 
von Neumann stable. 

5. The method is applicable to rectangular domains. Although it does not 
necessarily preserve symmetry along diagonals of rectangular domains, the 
procedure can be corrected if so desired (Ref. 20). 

A major attraction of the Douglas-Gunn Bcheme is that the intermediate 

* * * n +l 

solutions $ and $ are consistent approximations to <J> . Furthermore, 

D A kk n +l 

for steady solutions, $ = $ ** $ = $> independent of At. Thus, physical 

boundary conditions for can be used in the intermediate steps without a 

serious loss in accuracy and with no loss for steady solutions. In this 
respect, the Douglas-Gunn scheme appears to have an advantage over locally one- 
dimensional (LOD) or "splitting" schemes, and other schemes whose intermediate 
steps do not satisfy the consistency condition. The lack of consistency in 
the intermediate steps complicates the treatment of boundary conditions and, 
according to Yanenko (Ref. 14, p. 33), does not permit the use of asymptoti- 
cally large time steps. 

It is worth noting that the operator ^ can be split into any number of 
components which need not be associated with a particular coordinate direction. 
As pointed out by Douglas and Gunn (Ref. 15), the criterion for identifying 
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sub-operators is that the associated matrices be "easily solved" (i.e., narrou- 
banded) . Thus, mixed derivatives can be treated implicitly within such a 
framework, although this would increase the number of intermediate steps 
and thereby complicate the solution procedure. Finally, only minor changes 
are introduced if, in the foregoing development of the numerical method, 2A, 
and S are functions of the spatial coordinates and time, as well as $. 

Alternative techniques to the consistently split linearized block implicit 
scheme presented here have been proposed as general algorithms. We discuss 
two such methods, a predictor-corrector scheme due to Rubin (which we will 
refer to as P/C) and the Hopscotch algorithm of Gourlay and his co— workers. 
Although these authors have had success with these methods, it will be shown 
that in order to meet the requirements of the problem under consideration, 
they are not as versatile or as efficient as the consistently split linearized 
block implicit scheme previously discussed. 

Predictor-Corrector Method 

A predictor-corrector method has been successfully employed by Rubin and 
Lin for three-dimensional viscous flows in which diffusion is important in 
two directions (Refs. 5 and 21). Their objective in developing a compromise 
between explicit techniques-and implicit methods (ADI) was to eliminate viscous 
stability restrictions and to minimize CFL stability limitations common to 
explicit methods and to reduce the total work per time step by eliminating one 
of the block tridiagonal inversions required in the usual ADI procedure. In 
addition they desired a "symmetric" method that would easily treat symmetry 
conditions and other derivative boundary conditions. 

A comparison of ADI with Rubin's predictor-corrector method can be 
obtained by considering the model two-dimensional linear parabolic equation. 

u t “ (au yy + eu y ) + ( b u 2r + d u 2 ) (6) 

The P/C method reduces to the following expression 

{i - /3AtD y +2\ z /3b}[u**-u n ] ■ At{D y +0 2 }u n 
+ /3AtD z [u*-u n ] + 2 X. z /3b(u* -u n ) 
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where u denotes u at a previous iteration, u denotes u at the latest 

2 

iteration. D “ au + cu D ■ bu + du„ and X ■ At/Az . As a first 
y yy y z zz z z 

guess one of the following extrapolations is used; 

* D 

linear replacement u » u 

O(At^) Taylor scries u* - 2u n - u n “* 

The use of linear replacement achieves consistency only with many iterations 
or small time steps. 

The sequence of steps for solving Eq. (7) by the P/C method are: 

1. In Eq. (7) compute At(D y + D z )u n . 

2. Compute 8AtD z (u* - u n ) and 2X z 8b(u* - u n ) using u* from a previous 
iteration. 

3. Compute the coefficient matrix of (u** - u n ) . 

A. Obtain (u** - u n ) (requires one tridiagonal sweep). 

5. Repeat steps 2-4 until convergence, usually two to three iterations. 

Note that even for linear problems, iteration is necessary to obtain the 
desired accuracy. Thus at a minimum two to three tridiagonal sweeps and two or 
three explicit evaluations of the terms in step two are required. 

By comparison the Douglas-Gunn ADI procedure gives for the combined two 

steps 

[ I - /3AtD y ][V* - u n ] - A»[o y + D z ]u n + /3AtD z [u* - u n ] (8) 

* ** 
where u corresponds to an intermediate solution and u to the solution at the 

(n+l)st time step. Note that Eq. (8) differs from Eq. (7) in appearance only 

in' the underlined terms in the former equation. 

The sequence of steps in the Douglas-Gunn procedure for the two-dimensional 

problem is as follows: 

1. Compute the right-hand side to be used in the first sweep At[D + D_]u n . 

* 7 

2. Compute (u - u n ) , which requires one tridiagonal inversion, i.e., 

(I - 8AtD z )(u* - u n ) » At(D y + D z )u n . 

3. Compute (u** - u n ) which requires one tridiagonal inversion, i.e., 

(I - BAtD y )(u** - u n ) - (u* - u n ). 

4. Evaluate u n ^ <* (u** - u n ) + u n . 

The major effort is expended by both methods in evaluating the term 
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At(Dy + D z )u n . Neglecting the evaluation of the coefficient matrix of 
(u** - u n ) , two passes of P/C is equivalent to one double sweep of ADI. Hence 
by comparing the sequence of steps for both P/C and ADI we see there is no 
apparent advantage to the predictor-corrector method. Furthermore, due to the 
"explicit" nature of the P/C method in the " 2 " direction it possesses a CFL 
stability condition while ADI does not. 

If we were to consider higher order spatial approximations in both 
coordinate directions, each operator evaluation could involve a matrix inver- 
sion, (see section on spatial approximations). Hence the P/C method could be 
more costly than an ADI procedure. Finally the P/C method, to the authors' 
knowledge, has not been formulated for a full three-dimensional problem (with 
diffusion in all three directions, such as is required for the present problem) 
so that its applicability under such conditions is unknown. We see no 
advantage to the P/C method and therefore do not recommend it for the present 
problem. 


The Hopscotch Algorithm 


The so-called family of Hopscotch algorithms have been advocated by 
Gourlay and his coworkers (Refs. 22, 23 and 24) for the solution of multi- 
dimensional parabolic equations. Several variants of Hopscotch exist; for two- 
dimensional problems there is the fully explicit odd-even Hopscotch which 
resembles the DuFort-Frankel algorithm and the partially implicit line and ADI 
Hopscotch procedures which resemble the Peaceman-Rachford ADI method. These 
methods have lower operation counts than ADI schemes due to their partial 
explicit nature that, depending on the scheme, eliminates some or all of the 
matrix inversions. For instance, line Hopscotch requires only half the number 
of matrix inversions of a comparable ADI computation. However, the Hopscotch 
methods have stability restrictions and are only first order accurate in time. 

Numerical results of model linear scalar parabolic equations (Ref. 23) 
when the stability conditions are not violated, confirm the above conclusions, 
i.e.. Hopscotch is more efficient than ADI. However, recent results for the 
driven cavity problem (Ref. 25), which requires the solution of a coupled 
system of nonlinear equations, lead to contrary conclusions. The Poisson 
equation for the stream function was solved separately by a direct method 
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(Ref. 26), while the vorticity transport equation was solved by both ADI and 
odd-even Hopscotch. The ADI solutions were 20% faster than Hopscotch to 
obtain a converged steady state solution. Although Hopscotch requires fewer 
operations per time step than ADI, the CFL stability condition necessitates 
the use of a smaller time step which results in the noted Increased running 
times. 

More recently, Greenberg (Ref. 27) has added a number of new members to 
the Hopscotch family for three-dimensional parabolic problems. As with the 
two-dimensional procedures, these new members also have restrictive stability 
conditions. To date there has not been widespread use of Hopscotch-type 
schemes, in particular for coupled nonlinear parabolic equations, and their 
viability under such circumstances is still an open question. 

In general, Hopscotch owes its favorable characteristic to what Gourlay 
terms E-operators, of which the standard tridiagonal finite differences is a 
member. For other higher order spatial differencing, where the operators are 
tridiagonal but not E-operators, e.g., evaluations of spatial operators of the 
form au^ + bu x which even for an explicit temporal scheme involves a matrix 
inversion, the method loses most of its desirable features (see section on 
spatial approximations) . In addition the ability to handle coupled implicit 
boundary conditions is also not as flexible with Hopscotch. Finally, there 
have been claims that one of the attributes of Hopscotch is its ease of 
programming. While this may be true for model problems, for more complex 
problems, i.e., coupled systems of three-dimensional nonlinear parabolic 
equations, the basic logic for setting up the block inversions, whether in one 
direction as in Hopscotch, . or in three directions as in split LBI, is 
comparable. Hence, in this case, programming considerations should not greatly 
influence one's choice of method. 

The lack of versatility of the method, i.e., stability restriction, 
inability to incorporate higher order spatial methods, and results of more 
realistic problems (Ref. 25) lead us not to recommend the Hopscotch algorithm. 
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Spatial Difference Approximations 
Implicit Trldlagonal Finite Differences 

Q-R Operator Notation 

In this section implicit tridiagonal finite difference approximations to 
the first and second derivatives and to the spatial differential operator 
will be considered. The very versatile Q-R operator notation will be intro- 
duced, which allows as special cases a variety of schemes such as standard 
second order finite differences, first order upwind differences, fourth order 
operator compact implicit (OCI), fourth order generalized OCI and exponential 
type methods. Since all these schemes are of the same form, a single sub- 
routine which defines the difference weights is all that is required to 
identify the method, while leaving the basic structure of the program 
unaltered. Subsequently, the results of numerical experiments for a number 
of these schemes will be presented. 

The Q-R formulation allows for ADI methods and permits the treatment of 
systems of coupled equations, i.e., LBI methods. Although variable mesh 
schemes can be employed within the Q-R framework, it is believed preferable 
to use analytic transformations to obtain a uniform computational mesh, hence 
attention is restricted to uniform mesh formulations. 

The general concepts and notation for two point boundary value problems 
will be introduced and then the methodology entended to more general linear 
and nonlinear parabolic partial differential equations in one dimension. The 
extension to multidimensional problems will also be indicated. 

Consider the two point boundary value problem 

L(u) « a(x)u xx + b(x)u x + c(x)u - i (x) < 9 > 

with u(0) and u(l) prescribed. Derivative boundary conditions, although not 
treated here, can easily be incorporated into the framework of the Q-R 
operator notation. Let the domain be discretized so that Xj “ (j-l)h, 
j - i, 2,..., J + 1, and Uj % u(x.j), -v u x (x^), u^x^) and h = 1/J 

is the mesh width. The numbering convention was chosen here to be compatible 
with FORTRAN coding. 
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Without loss in generality for a(x) i 0, Eq. (9) can be divided by a(x) 
bo that we nay treat instead the following equation 
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L(u) " u xx + b (x)Uj. + c (x)u - f(x) 


( 10 ) 


where 


b(x)-b(x)/a(x), c(x)-c(x)/a(x) ond f(x) - T(x)/a(x) 

The spatial differential operator is identified as 

* u xx + b(x)u +c(x)u 


(ID 


Substituting the finite difference approximations to the first and second 
derivatives 


*2h' U j " 2h * F j * u x (x j ) + OCh 2 ) 


D + D - ■■ U j-r2U j+ U i+l 

h 2 V' * s j - u xx ( Xj ) + 0 (h 2 ) 


( 12 ) 


(13) 


into Eq. (10) and rearranging, we obtain 


L(u) ~ Sj + bj Fj 


u h ♦[«, - ipr] u j ♦ [ tJt * -If] u,., - 


[' - ] u i-i + [^cj-aJoj + [ i + -^-]u jM ■ h=fj 

where Rc^ = hb^ is the cell Reynolds number. 


( 14 ) 
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Equation (14) can be generalized by introducing operator format, i.e., 

r j u J-« + r ‘ u J + r J + U J4| • h2( f J-. + q* 1 fj + q} f ,♦,) (15) 

where the superscripts (-) minus, (c) center, and (+) plus indicate the 
difference weight that multiplies the variable evaluated at the (j-1), (j) and 
(j+1) grid points, respectively, and where the Tj's and q^'s for grid point j 
are functions of h, bj_ lf b^ , bj + ^, c j and c j+]/ Comparing Eqs. (14) and 

(15) we can identify the r^'s and q^'s, viz.. 


rj - I - RCj /2 
r j - h z Cj - 2 


qT * o 


fj ■ t + Rcj /2 


q| ■ o 


We now define the tridiagonal difference operators Q and R 

r [ u j] ■ vh + r 5 u j + r j + u J+ , 
o[ f j] ■q] f )., +qjfj + qJ^fj+i 

Noting that L(u) = f and substituting Eq. (17) into Eq. (15) we obtain 


R [Uj] = h z o[ l. ( U ) j ] = h z o[ f j] 


Alternatively by employing the inverse operator Q an expression for L(u)j 
can be obtained 

L(u), » -^-Q-'RUj 
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For standard central finite differences Q ■ Q~^ = I, the identity matrix, 
(the spatial operator is given explicitly in terms of U, U, and U..,) so 
that nothing was gained in obtaining Eq. (19). However, in general, for higher 
order methods Q is tridiagonal and Q - ^ is a full matrix. Hence Eq. (19) gives 
us a means of expressing the spatial operator for a wider class of difference 
approximations. The formalism in Eq. (19) is also applicable for first and 
second derivatives appearing alone (cf. Ref. 28). It must be pointed out however 
that Eq. (19) is not the most general formulation since the compact implicit 
formulas cannot be combined to yield a single scalar equation relating the 
spatial operator to the function values (Ref. 28). 

In the next section a method due to Berger et al., (Ref. 29) is described 
that enables one to construct fourth order tridiagonal methods with certain 
desirable properties, i.e., evaluate the q^ and r^ coefficients. 

Generalized Operator Compact Implicit Schemes 


Given 

I /ill el, 

'XX ■ " '“X 

an expression relating L(u) and u is sought in the form 


L(u) » + b(x)u„ + c(x)u 


I 


RUj ■ Q ( Lu)j + T j 


( 20 ) 


where x. is the truncation error and Q and R are tridiagonal displacement 

J 4 

operators. The maximum accuracy attainable is fourth order, i.e., x^ ■v 0(h ). 

Expanding Eq. (20) in terms of 9j’ C,+ anc * r j' C ' + we 


I 

If 


RUj- Q(Lu >i * [ r j u j-i + r j u J + r f u j*i] 

" [ q j < Lu)j_, + qj (Lu)j + qJ(Lu)j +I ] 


( 21 ) 
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A Taylor series expansion yields for t 


J 


ORIGINAL PASS I". 
OF POOR QUALITY 


Tj • T° u{ Xj) +TjU ,,, (Xj) + Tj 2 u (c, (Xj) + T j s u ,s, (*j) 
+ TjV 4) (Xj) + Tj°u <5, { Xj) + T*u , 8 , (Xj) + 0 (h°) 


where superscripts in parentheses denote derivatives with respect to x, and 
where 

Tf - “^[( 1-7 +r, C + r]) + h 2 (q,^ J . l +q*c J + q}c j4 ,)] 

T j' * 4 ” [ (r f - rp-h(qj-bj., + qjt> J +qJ’b J+ ,)-h 2 (q]c Jt| -qJc,. 1 )] 


Tj 2 - Y (r f + r j" } - ( q j” + qf + qp 

jz 


(23) 


- h(qj*bj +1 - q[b H ) - -7- (q*c j+1 + qjcj-.) 


v , V -2 

Tj • n 


if * ‘-""'P ' + C-D—qj-bj.i 1 


[ q j* + (-1)%-) + ti 2 (q,*C|„ + (-l)'qjc,.,) 


v * 3, 4, 5,6 


For second order central finite differences we set T° - T 1 - T 2 - 0. This 


yields, when q^ - 1 and q~ « q* » 0, the following relations 


c 

r i ■ 


-(rj + rj*) + h Z ( q j"Cj.| + qjCj + qj + Cj„) 


Tj* - Tj” ■ hbj ■ Rcj 

rf +r P 2 
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which recovers Eq. (16), i.e.. 


original '■‘J 

OF POOR Q'J.VJlV 


r] ■ -2 + h z c 


r“ - I - Rcj /2 
rj* ■ I + RCj /2 


To obtain the fourth order operator compact Implicit scheme ve again set 
X® » " 0 to obtain throe expressions for rj ,c,+ in terms of the 

qj ,c,+ , (note that q~, q + 4 0 and q c is not necessarily unity), i.e.. 


r 5 ’ * (r j + + r r } + h2(q f c j-i + q j e<: i +q j* c j + . ) 


(24a) 


r j 4 - Fj" * + RcjqJ + Rcj + 1 qj + + h a (q*Cj +| + qfcj.,) (24b) 

r* + r“ • 2(q‘ +qj +qp +2[Rc J+) q^ - Rc^qj"] + h z (q*c H + q^Cj.,) (24c) 


3 4 - c + 

Now T and T must be set to evaluate q ’ ’ . The standard Swartz OCI 

3 4 ^ 

method requires T “ T “0, 


_ (r *-r-)--i-[Rc j , | q j + + Rc^.q"] - [qf-qf] -0 
-h 2 (qp J+| -qj-c H ) 


(24d) 


I 

IT 


(r * +r) ‘ IT [ Rc jW ' Rc j-, q j] “ T [ q J* +q J _ ] "° 

“ h 2 (qj +c j+ , + qj'Cj.,) (2«e) 
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and results in a leading truncation error term of (Au^ + Bu^)h\ Substi- 
tuting (24b) and (24c) into (24d) and (24c) r and r can be eliminated and a 

- c + J J 

system of two equations in q^ and q^ with qj as a parameter is obtained. The 
parameter q^ is proportional to the determinant of the system. The values of 
q * * and r* ’ are presented in Table I. 


’j 




As shown in Ref. 28, a cell Reynolds number stability condition exists 
for the Swartz OCI scheme, i.e., for Rc £ /II nonrealistic or oscillatory 
solutions will be obtained. In order to eliminate this restriction one can 

3 4 (3) (4) 

relax the conditions T ■ T » 0, and allow the coefficients of u and u 

to be of O(h^) . 

By expanding q 


•»c,+ 


j 


in a series in Rc 


. — .c ,♦ 


I X m' C,+ Rc 

m*o 


m 


(25) 


.. c *f 

12 parameters, X m ’ ’ a ■ 0, 1, 2, 3, are introduced. The equations for 
T^ » O(h^) and T^ ■ O(h^) yield 5 linear relations, leaving for disposal 6 
"free" parameters plus a factor. 

These parameters can be set according to some criteria that would yield 
certain desirable properties for the difference equations. The following 
constraints are prescribed 


Qj* >0, qj" >0, qj C >0 

q ] Cb | * q f b ]-| +qf b j»i 


(26) 


and h is sufficiently small so that 

I0bj-bj.,-bj +1 > 0 and hcj +1 /bj +| <2 
for J-2, • ,J and Cj£ O 
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These conditions assure that R Is diagonally dominant and Q is invertible for 
all R fi . Further details arc given in Ref. 29. The significance of this 
approach is that one can construct a scheme possessing certain desired proper- 
ties employing a set of preassigned rules. This is contrary to usual practice, 
in which a scheme is chosen, and then its properties are determined. It is 
important to note that the derivation of the q and r coefficients is not a 
trivial task. In computational effort it is also not cheap. This point will 
be discussed in greater detail in a later section. 

The q and r coefficients for the generalized OCI scheme described in 
Ref. 29 are given in Table II. In this report comparisons are made with another 
generalized OCI scheme, whose coefficients are given in Table III. Numerical 
experiments indicate that these two schemes are comparable, differing only in 
the magnitude of the truncation error. 

Exponential Type Schemes 

Another family of schemes that can be expressed in Q-R operator notation 
are the so-called exponential methods. The idea, originally due to Allen 
(Ref. 30) (independently derived by Il'in (Ref. 31) and McDonald (Ref. 32)) 
and employed by Dennis (Ref. 33) , is to set the difference weights so that the 
numerical solution is equated to the analytic solution for the locally frozen 
constant coefficient equation. Allen (Ref. 30) and Il'in (Ref. 31) considered 
the homogeneous constant coefficient equation, 

L(u) • U xx + bu x • 0 (27) 

so that the difference approximation was set identically equal to the analytic 
solution. 

The analytic solution of Eq. (27) has the form 

. . -br 
U * A + Be 

where A and B are determined from the boundary conditions. 

The Q-R operator formulation of Eq. (27), 

r f U H * 'l‘ U i * '!%> ■ 0 (26) 
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possesses a fundamental solution of the form y^. Solving for y and employing 
Eq. (24a) we obtain 

H- r j* /r f 


With the aid of Eq. (24b) and setting y^ to the analytic solution e -b J x we 
obtain 


ft) . e' b J x - e' b i hJ • e* Rc J - r.Vrf 


I ' * J 

With the condition r* - r” - RCj (T* - 0 in Eq. 23), we can define r” ,C,+ 


as 


r" • Rc e* Rc J / (I - e* Rc J) 

r* - Rcj /(I -e* Rc J) ( 29 ) 


» -f C 

and q “ q ”0 and q - 1 where we have allowed b^ to vary. Equation (29) can 
be rearranged to yield an alternate form 

L(u) j ■ ~r~ t co,h “iH D *°- u j + b j D o u j (30) 

where D + , D_, D q are the forward, backward and central first difference oper- 
ators, respectively. 

This method is second order accurate for Rc 'V/ 0(1) and becomes first order 
accurate as Rc + ■ where the scheme reverts to first order upwind differencing 
(in Eq. (28) rj -*■ 0, r j -* Rc^ , rj -Rc^ for b^ > 0 to give Rc j ( U j+i “ U j) " 
Another exponential scheme which is uniformly second order accurate was 
developed by El-Mistikawy and Werle (Refs. 34 and 35). The "exponential box 
scheme" which is incorporated in their solution of the boundary layer equations 
with strong blowing, is based on a spatial operator of the form given in 
Eq. (10). Berger et al., (Ref. 36) derived the counterpart for an operator of 
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the form given in Eq. (27). The q and r coefficients are presented in 
Table IV. Although this scheme reverts to second order upwind differences as 
Rc -*■ <•, it does not possess a maximum principle analogous to the ordinary 
differential equation it is approximating as does the exponential scheme of 
Allen (Ref. 30). 

The Allen exponential scheme can be applied to differential equations of 
the type given in Eq. (10). Substituting the Q-R representation for L(u), we 
obtain 

■Ijr q -i r[uj] + CjUj ■ fj < 31 

Multiplying through by Q, and combining terms, the difference approximation to 
Eq. (10) becomes 


r[ Uj] + h Z Q [ Cj Uj] - h Z Qfj 

or 

(R +h Z QC J )[u j ] • h 2 0fj (32) 


Note that the only difference between Eq. (32) and Eq. (18) (where Cj 5 0) i:; 
the coefficient matrix multiplying Uj. Hence the methodology is unaltered. 

Application of Q-R Operator Schemes to 1-D Parabolic Equations 
Consider the one-dimensional linear parabolic equation 


u, * au w + bu x + cu + d 


(33) 


with appropriate boundary conditions and initial conditions, where a, b, c, d 
could be functions of x and t. Dividing by a i 0, identifying L(u) ■ u^ + bu x , 


If Eq. (33) were nonlinear, assume it was linearized by the method described 
in Appendix B. 
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and employing a two point temporal difference scheme, we obtain 



L J n * 3 (fl U| n "+(l-/3) U| ") + d J n * /3 


where 8 ” 1 reduces to a backward difference and 8 “ 1/2 to the Crank-Nicolson 
approximation. 

Letting 


and with some algebra we obtain. 



2 

where X - A t/h . 

Again, Eq. (34) is general and permits second order finite differences, 
upwind differencing, exponential type schemes, and OCI schemes. Aside from the 
evaluation of the Q and R operators, the problem is no more complicated than 
standard second order finite differences. Results of model problems can be 
found in Ref. 28. 


Application to Coupled Nonlinear Parabolic Equations 

Given a system of m nonlinear parabolic equations in m unknowns, 


m 

X 

1=1 


. n+i n, 

(Ujj -Ujj) ..n+/ 3 , 

n +/3 At " N i H (u 1 ,u 2 ,....u m ,x | ,x 2 ,x 3 ,t) 


'!] 


j - 1,2,..., J+l 
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where n" P is a nonlinear nonconservative spatial operator, the Q-R formalism 
carries directly over provided that for any equation only one independent 
variable is operated upon by the differential operator. For example. 


I 



is allowed since x derivatives of u only appear, while 

I 

alu.'wTv) U ‘ * U « + b(u ’ v ’ w)u * + c(u > v » w) + d(u,v,w)w x 

is not allowed since x derivatives of both u and w appear. The modified 
unsteady Navier-Stokes equations for the three-dimensional time-dependent 
boundary layer, when written in quasllinear form, fall within the class of 
allowable differential operators. Thus for the problem being addressed 
in the present study the OCI schemes are applicable. 

Multidimensional problems and/or more general equation forms can usually 
be accommodated by a splitting procedure, which reduces the differential 
operator to a sequence of one-dimensional problems which have the appropriate 
allowable form. However, as with standard finite differences, special proce- 
dures must be applied to cross derivative terms, e.g., extrapolation or lagging 
at the previous time step or increasing the number of intermediate steps in the 
splitting. 

In Appendix B an example of the above procedure is presented for a 
coupled system of one-dimensional nonlinear parabolic equations possessing 
"boundary layer-like reverse flow" behavior. 

Basis Functions 

A convenient approach for developing numerical procedures for the solution 
of partial differential equations is founded on a basis function representation 
of the dependent variable, i.e., in one dimension, 

u(x) * Ia,B:(x) 

I 1 1 
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where u(x) is represented by a linear combination of suitably chosen basis 
functions, B^(x), and the a are determined by some set of constraint relation- 
ships. The advantage of the basis function representation is that once the 
are determined, one can obtain at very little additional cost the function 
and its derivatives anywhere in the domain. 

The major concern here will be with the B-spline basis functions. However 
by way of comparison some of the characteristics of piecewise Lagrange poly- 
nomials and piecewise Hermite polynomials will also be given. In the following 
basis functions are sought, polynomials of order k (degree k-1) , that possess 
certain smoothness properties, and the computational efficiency of these 
functions for the solution of differential equations are investigated. 

The simplest functions are the piecewise Lagrange polynomials, which can 
be computed by the cardinal basis functions. Consider a grid numbered from 
j ■ 1 to J + 1 (to be consistent with FORTRAN coding) so that one is consid- 
ering J + 1 grid points and J intervals. Over each interval consisting of k 
knots the basis functions are polynomials of degree k-1 which are equal to 
unity at one particular knot and zero at the other knots. The dimension of 
the piecewise Lagrange polynomials (the number of basis functions) is J + 1 
(independent of the degree of the polynomial) and thus J + 1 constraints 
must be satisfied to determine all the a^. However, in order to evaluate a 
function and its derivatives at a particular location, only k basis functions 
distributed over k adjacent grid points are required. In addition Lagrange 
polynomials are only C° at the end points of each sub-interval, again 
independent of the degree of the polynomial and thus allows for jumps in the 
first derivative there. Due to the purely interpolatory character of piece- 
wise Lagrange polynomials the are nothing more than the functions evaluated 
at the appropriate knots. 

The interpolation problem reduces to an explicit procedure ns does the 
evaluation of the derivatives at the internal knots. When incorporated with 
the method of collocation (the analog of interpolation for the solution of 
ordinary differential equations) , the piecewise Lagrange basis functions used 
locally on a uniform mesh (with k odd) recovers the standard centered finite 
difference approximations. Thus a quadratic polynomial leads to a three point 
formula and the inversion of a tridiagonal matrix of order J + 1 and a quartic 
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polynomial leads to a five point formula and the inversion of a pentadiagonal 
matrix of order J + 1. 

In order to obtain smoother representations additional constraints must be 
applied on the derivatives of the approximating functions. This leads to 
more complicated basis functions. For example, the family of Hermlte poly- 
nomials requires that for a (2m - 1) degree polynomial (m i 2) on two adjacent 
knots, m interpolatory constraints must be satisfied for the function and its 
(m-1) derivatives. The approximation is still local, i.e., given a function 
and an appropriate number of derivatives on two adjacent knots, the function 
and its derivatives can be computed implicitly in the interior of that domain. 
However, the solution of a differential equation by the method of collocation 
involves the inversion of a matrix of order m(J + 1) of bandwidth (3m - 1) 

(■ n + m where n «• degree of polynomial) . Thus in the process of obtaining a 
C n ^ representation a substantial increase in labor in comparison to Lagrange 
polynomials has been incurred. Further details can be found in Ref. 37. 

B-Spllnes 

Another family of functions are the polynomial splines, i.e., polynomials 
k-2 

of degree k-1 that are C . One would expect that the additional smoothness 
of the spline approximation would translate into a better behaved solution - 
perhaps even more accurate - for “smooth" functions, than a Lagrange or 
Hermite polynomial of the same degree. Numerical experiments have shown this 
to be the case in many instances. However, the additional smoothness 
constraints that must be specified not only change the nature of the approxi- 
mation, but may add considerably to the computational effort, depending on the 
manner in which they are applied. 

There are several approaches for representing polynomial splines. The 

first method involves specifying the polynomial and the smoothness conditions 

£ 

separately, and then solving the imposed constraint equations simultaneously 
for the undetermined coefficients. This approach has successfully been 


The constraint equations may either be interpolatory conditions or a 
differential equation. 
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employed by Kendall and Bartlett (Ref. 38) In their solution of the chemically 
reacting boundary layer equations, and Murphy et al., (Ref. 39) for the solu- 
tion of the boundary layer equations. With the aid of a skillful partitioning 
of the resulting banded block matrix, the above authors were able to solve 
for the function values directly. However, by so doing they lost ready access 
to the interpolatory polynomial representation, which is one of the acknow- 
ledged benefits of a spline technique. An alternate approach described in 
Ahlberg, Nilson and Walsh (Ref. 40) and adapted by Rubin and Khosla (Refs. 41 
and 42) in their solution of the boundary layer equations also employs a 
polynomial formulation. However, the additional continuity and smoothness 
conditions are used to transform their equations into expressions relating the 
function values and their derivatives. Their "spline equations" are solved 
simultaneously with the appropriate constraint equations, which entails for a 
single scalar differential equation the inversion of 2 x 2 block tridiagonal 
matrix to obtain fourth order accuracy and a 3 x 3 block tridiagonal matrix 
to obtain sixth order accuracy. 

In order to avoid inverting block banded matrices that arise even in 
the solution of scalar differential equations a third alternative is 
considered, the B-spline representation, which has the appropriate smoothness 
built into the functions themselves. For the solution of the coupled Navier- 
Stokes equations this procedure has a definite advantage in that the order of 
the block submatrices, typically three or four, can be reduced by one. A 
detailed discussion of B-splines, their construction ar.d their mathematical 
properties are given in Refs. 37 and 43 and will not be repeated here. 

Instead, the properties of .B-splines that make them attractive for the solution 
of partial differential equations will be emphasized. In the following 
discussion, without loss in generality, a uniform mesh will be considered. 

The normalized B-spline (the sum of the basis functions at any knot is 
unity) of order k (degree k-1) is bell-shaped and spans k + 1 knots. Within 
the interval the B-spline is positive while outside the interval. Including 
the end points, it is zero. For knots of multiplicity unity the first k - 2 
derivatives of the B-spline are zero with a jump possible in the (k - l)st 
derivative. In general, for a knot of multiplicity k - v the vth derivative 
of the B-spline is discontinuous at the point with all lower order derivatives 
continuous. Due to the compact support property of the B-splines at most k 


35 





B-splines are nonzero over any Interval, while at a knot k - 1 B-splines are 
nonzero. This property enables one to reduce the bandwidth of the matrices 
that are obtained in the solution of differential equations, while retaining 
the scalar structure. 

DeBoor (Ref. 44) has shown that for B-splines of order k the dimension 
of the basis is J + 1 + k - 2, where it is assumed that the knots are of 
multiplicity unity and there are J + 1 nodal points (J intervals) in the 
domain. Thus for cubics (k ■ 4) there are J + 3 independent basis functions 
and for quintics (k “ 6) there are J + 5 basis functions. Hence J + 1 + k - 2 
conditions are required to fully specify the B-spline representation. In view 
of the above, for a B-spline of order k, in order to determine all the a^, the 
resulting matrix is scalar, is of order J + 1 + k - 2 and has a bandwidth of k. 
After the are determined, the evaluation of the function and its k - 2 
derivatives each require at most k additional multiplications. 

For the pure interpolation problem 



K - 3 



j ‘ I J+l 

G - I, .... K-3 


where superscript i, designates a derivative of order l, and k even, Prenter 
(Ref. 37) shows that a unique function exists, the polynomial spline of order 
k. What is apparent from the statement of this problem is that the B-spline 
interpolating polynomial is global in character, being "tied together" by 
the derivative constraints at the end points. In contrast to piecewise 
Lagrange and Hermite polynomials where coefficients were determined for each 
subinterval by the appropriate interpolatory constraints, the evaluation of 
B-splines requires the inversion of a matrix. In addition the a^'s that are 
determined are not equal to the values of the functions and its derivatives 
at the knots as was the case for the Lagrange and Hermite representations. 

It Is in part due to this global character that B-splines derive their 
smoothness properties. 
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DeBoor (Ref. 44) has presented a computer program to determine the 
normalized B-splines of order k < 20, which allows for nonuniform meshes and 
for discontinuities in lower order derivatives at the knots (knots of multi- 
plicity > 1) . He has also conditioned the B-splines near the boundaries to 
simplify the setting of boundary conditions. Alternatively the procedure 
suggested by Prenter (Ref. 34) could be employed (for a uniform mesh with 
simple knots) , but that would require additional algebraic manipulation to 
modify the coefficient matrix so that one can obtain the appropriate band 
structure. Since both methods are equivalent (the uniqueness of the B-spline 
representation) the choice is a matter of convenience. Prenter's approach 
employs the same basis functions throughout the domain, which reduces the 
memory and computer logic requirements but increases the algebraic manipula- 
tion, in particular for higher order splines. The more general DeBoor 
approach requires more computer memory, but reduces the preprocessing required 
of the user. In the subsequent work DeBoor' s approach was used. 

The accuracy of the B-spline representation is now considered. For a 
function u that is sufficiently smooth, with derivatives of order k that are 

k 

continuous, u e C (a,b), the B-spline representation of u , s (of order k, 
k even) , Prenter gives the following error estimates 

llu-sll*, - 0(h k ) 
llu'-s'll„ - 0( h k "') 
llu"-s"ll 0J - 0(h k " 2 ) 

Hence, to solve a second order ODE by collocation using cubic splines (k «* 4) 

2 

we can expect 0(h ) accuracy but require only a simple tridiagonal matrix 
elimination. However, to obtain greater accuracy larger bandwidth matrices 
are required. Consideration of this point is given in the following section. 

Application of B-Splines to the Solution of Differential Equations 

Since the B-spline representation will be employed as a spatial approxi- 
mation, its properties, order of accuracy, efficiency and spatial stability 




behavior can be investigated by considering its application to the solution 
of two point boundary value problems. 

In this section B-spline basis functions of order k are applied to 
the solution of ordinary differential equations of the form 


L(u) - U xx + bu x + cu « f 


-[0.1] 


with boundary conditions 


, U(l) 


Several different approaches can be employed for the solution of Eq. (35), 
depending on the degree of accuracy desired and the level of complexity one is 
willing to accept. The techniques of collocation, Galerkin, subdomain, 
least squares, etc., can be viewed as special cases of the method of weighted 
residuals (of Crandall (Ref. 45) for instance). More recently Murphy (Ref. 46) 
has shown the relationship of these methods to orthogonalization processes 
and has characterized them as generalized Galerkin techniques. Here we will 
review the relationship of the different approaches to the method of weighted 
residuals and indicate how- they translate into computational effort. 

The basic idea of the method of weighted residuals is to choose a trial 
function for the independent variable 

u(x) " ZajB kj (x) (36) 

where are the basis functions (here B-splines of order k) and are 
coefficients to be determined. The derivatives that are given by 


u'(x) - I ctj B kj (x) 


u"(x) ■ Z a ( B k| (x) 


are substituted into Eq. (35) to yield 
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L(U) 



L(u) ■ L^a, B ki (x)» ^ajLB k |(x) 
B k j(x) + b( x) B^tx) +c(x)B k |(x 


>} 


We now require that the residual error in some Gcnse vanishes over the domain, 
i.e.. 


X / [lu(x) - f(x)] V/ m (x)dx « 0 


(39) 


where is a suitably chosen weighting function. The simplest method is to 
force the residual to vanish at each of the nodal points in the domain. This 
method of collocation is analogous to the interpolation problem which requires 
the approximating function and/or its derivatives to be equal to the true 
function at the nodal points. 

In this case W is the Dirac delta function and we recover the system 
□ 

of equations 


Lu(x m )- f{x m ) - 0 m ■ |,... , J + | 
or 

^ a i[ B ki U i ) + b(X j )0 ki (X ) ) + C(x J )0 KI (x j ) ] ’ f j + l (40) 

For B-splines of order k there are J + 1 + k - 2 basis functions 
(dimension of the B-spline basis) so that in addition to the collocation 
relationships at each of the nodal points, including the boundary points, k - 2 
supplemental conditions arc required. For cubic B-spline (k ■ 4) the two 
boundary conditions are sufficient to close the system. Hence a linear system 
of equations of order J + 3 and bandwidth k - 1 - 3 (tridiagonal) is obtained. 
Note that the problem is no mere difficult than the standard second order 
finite difference case, except that J + 3 equations are solved Instead of at 
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k-2 

most J + 1 equations. The overall accuracy of the method is 0(h ) which 

corresponds to the accuracy of spline approximation of the second derivative. 

Once the a. 's ore determined, u, u* and u" can be obtained from Eq. 07). 

* I 

For higher order splines, e.g., k ■ 6 (0(hj accuracy) an additional k-2"4 
conditions are required to close the system, which reduces to specifying two 
aore conditions than are available, i.e., the two boundary conditions. If one 
has information concerning the behavior of the differential equation at the 
boundaries (values of derivatives there), then that can be used. However, 
in general this is not always the case, so it was decided to collocate at two 
points that were noncoincidcnt with the knots. These locations were set at 
X ■ h/2 and at x ■ 1 - h/2 solely for ease of computation (modifying the 
matrix), but insight into the behavior of the differential equation could be 
used to choose the collocation points, i.e., regions of rapid change. It is 
important to note that these collocation points are located in a subregion 
between two knots. There is no need to add grid points (knots). 

The bandwidth of the resulting matrix has now been increased to 5 in 
order to obtain an O(h^) solution. If the coefficients, b, c and f are not 
known at x - h/2 and x ■ h/2 then these values must be obtained by interpola- 
tion. 

It can be shown that more accurate results can be obtained (one order 
better) if derivatives are introduced as unknowns, i.e.. 


u " fa, 8 ,,, 

u' - I/3jB ki 


(41) 


A system of two equations in 


Qj and 8^ are thus obtained 


L(u) r I {/3 i [Bki (x J )+b i B kl (X J ) ] +a i C J B kl (x ]>} ‘ f J 


and 


?fe B kl (X j ) " Q i B kf X j ) }‘ 0 


(42) 


(43) 
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However, a 2 x 2 block matrix 1 b now required to be Inverted and boundary 
conditions on u' must also be supplied. 

Calerkln Method 

If we choose the test function W n to be from the same family as the trial 
function, then the weighted residual becomes analogous to an orthogonality 
condition. Tills leads to the following relationship, 

[ Z a ! (B Hl +b(x , )B kl +c( *) )B kl>-f,Km dx - 0 (44> 

and we obtain (J + 1 + k - 2) equations which is exactly equal to the dimension 
of the basis. However, the bandwidth of the coefficient matrix for the a^'s 
is increased to 2k - 1 (7 for cublcs and 11 for quintics). This result can be 
obtained by focusing our attention at the B-6pline test function centered at 
knot j - m. It will contribute to the integral in the interval (m - k/2, 
m + k/2) where it is nonzero. However, the B-splfncs centered in the interval 
lm - (k - 1) , ra + (k - 1) ] will also contribute to the integral. Hence the 
B^’s will span (k - 1) + (k - 1) + 1 » 2k - 1 knots and will contain 2k - 1 
entries. 

The Galcrkin scheme can be shown to be O(h^) accurate for cubic B-splines. 
However, as compared to the method of collocation the bandwidth has increased 
from k - 1 to 2k - 1 so that the procedures would be much more costly. Further- 
more, the Galerkin scheme requires four nontrivial integrals which necessitates 
an integration scheme consistent with the order of the method, e.g., for cublcs 
a fourth order Simpson's rule. Since the B-splincs are not orthogonal in the 
sense that integrals of products of basis functions are not zero, no simpli- 
fications exist. 

Since in general one is not solving a problem that stems from a variational 
formulation, one would be Justified in searching for test functions from a 
different family which could simplify the integrations and reduce the bandwidth 
of the resulting matrix. One Buch technique, the subdomain method (Ref. 45) 
or what Murphy (Ref. 46) terms the generalized Galerkin procedure, employs a 
unit step function as the test function. This formulation yields the following 
system of equations 
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Since only J equations arc derived for the (J + 1 + k - 2)a^a, k - 3 additional 
relationships arc required. For cubics, thiB leaves, aside from the tvo 
boundary conditions, one as yet unspecified condition. Murphy (Ref. 37) and 
Bartlett and Kendall (Ref. 38) in their viscous flow solutions specified an 
additional derivative boundary condition at the edge of the viscous layer. 
Although their condition was exact mathematically for the problem they 
considered, for more general equations such conditions could overconstrain 
the solution. Hence it is felt that collocating at some point, i.e., one of 
the boundaries, would be helpful and would not deteriorate the order of 
accuracy of the solution (this point is further discussed in the section on 
numerical results). 

The bandwidth of the resulting matrix has now been reduced to k as 
compared to 2k - 1 for the standard Galerkin scheme while still retaining 
0(h ) accuracy. However, the integration scheme (now over two adjacent grid 
points) needs to be appropriately handled in order to achieve the desired 
accuracy of the method. In order to employ Simpson's rule for cubic B-splines, 
values of the coefficients b, c and f are required at the intermediate points 
Xj + ^ 12 ' If they are not known analytically, interpolations would be 
required which would increase the total computation time. 

The advantages of the cubic B-spline generalized Galerkin procedure are 
its decreased bmdwidth of four and its applicability to treat equations in 
conservative form. However, it has the disadvantage of being sensitive to 
the type of boundary condition used to satisfy the extraneous condition 
(cf., section on numerical results). 

The advantages of the quintic B-spline collocation procedure are that 
boundary conditions are easy to apply and it possesses the lowest truncation 
error for a given grid spacing. However, it is less efficient than the 
generalized Galerkin method due to its larger bandwidth of five. At the 
expense of increasing computer memory requirements, quintic B-spline collocation 
can be made more efficient. 

Both methods have the advantage over nonbasis function 3cheme3 of allowing 
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one to easily obtain derivatives of the solution vector as well ns the solution 
at any point in the donain. Their main shortcoming is in the cell Reynolds 
number stability condition (cf., section on numerical results) which manifests 
itself as oscillations in the solution. In view of the above a blanket 
recommendation cannot be given. Only after an analysis of the problem under 
consideration is made and the goals of the solution (accuracy and efficiency) 
have been stated can one recommend either of the schemes. 

Time was not available to program either Rubin's method or the Murphy- 
Kendall-Bartlett technique and thus a rigorous comparison with the B-spline 
formulation presented here cannot be given. However, it would appear that 
both B-spline procedures are competitive with these other spline variants, 
and they do have the advantage of reducing the order of the blocks for coupled 
systems of equations. 

A comparison of the B-spline algorithms with the Q-R operator techniques 
is given in the section on numerical experiments. 

Application of B-spllnes to 1-D Parabolic Equations 

* 

Consider the one-dimensional linear parabolic equation 

u, - au xx + bu x + cu + d ( 46 ) 


with appropriate boundary conditions and initial conditions, where a, b, c, d 
could be functions of x and t. Employing a two point temporal difference 
scheme, we obtain, 



0 ♦ (i-0 )„;■"> ♦ b J n * e (/3u J '"*' ♦ d-0) U ; n ) 

* c"‘ 0 <0u + (l-/3)u, n > + dj"** 3 


If Eq. (46) were nonlinear, assume it was linearized by the method described 
in the appendix. 


43 



where 0 - 1 reduces to a backward difference and 0 ■ 1/2 to the Crank-Nicolson 
approximation. 

The B-cpline spatial approximation can be combined directly with a 
temporal discretization scheme to solve Eq. (46). Here collocation and the 
generalized Galerkin procedures are considered. Although no numerical 
experiments were run with these methods, they are presented here to indicate 
how they can be applied. 

1. Collocation 

Substituting Eqs. (36) and (37) into (47) and rearranging one obtains 

? °r{w • + »rv«,> W8] 

+ C j n+ \i (x i ) ]}" A,d j n+/3 ’ RHS " j-',2 J+l 


where 


RHS 


4I(I-S>[ 


+ b ,"*' 2 3 


/ n , . n+/3 n l 
u ) + C J u j J 


+ u, 


is known from the previous time step and 

a ,"*' 3 . So,"*' + (I-S)o," 

and similarly for 


n+/3 

Dj , Cj 


and d"*^ 


2. Generalized Galerkin (cubic B-splines) 

Integrating Eq. (48) over (j, j + 1) one obtains 
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/ Xi *'[ K*' k.(x,) - * b, n Ai'<*,> 

X J 

+ c J n+/ 3 a kl (x J ))} - Atd^jdx » / Xj+ ’ (RHS n )j dx 

X J 

By the use of Simpson's rule the integrals can be evaluated in terms of the 

function values at x., x and x If the values at the "half points", 

n+B ,n+B n+B J j (n+B d ‘ . . ... ,, . . 

a j+l/2' j+1/2’ C j+l/2 and d j+l/2 are n0t known analytically they must be 

obtained by an interpolation routine. 
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RESULTS 


Numerical Experiments 

Spatial schemes can be analyzed by considering two point boundary value 
ordinary differential equations. The properties of these schemes, accuracy 
versus grid spacing, running time for a given accuracy, ease of programming, 
storage requirements and cell Reynolds number effects can easily be obtained 
without the encumbrance of temporal discretizations. 

Here the following differential equation is considered 


U *X + bU X " f 


(49) 


with appropriate boundary conditions specified. 

The first type of problem investigated was the homogeneous constant 
coefficient case 



u xx + bu x ■ 0 

u(O) • 0 , u( I ) ■ I 


(50) 


which has as its exact solution 


U exact 



(51) 


The aim of these numerical experiments was to verify the convergence rates 
of the various methods and the effect of cell Reynolds number on obtaining 
physically meaningful solutions. The convergence rates were verified by 
considering the case of b *» 1 for 10, 20 and 40 intervals, so that effects of 
cell Reynolds number, which ranged from .025 to .100 would be insignificant. 

In Table V these results are presented for the generalized OCI schemes (0CI-G1, 
(cf.. Table III), 0CI-G2, (cf.. Table II)), cubic B-spline collocation, (CBS-COL), 
quintic B-spline collocation, (QBS-C0L) , and cubic B-spline generalized 
Galerkin, (CBS-GAL) . All the schemes except for CBS-COL have fourth order 
convergence rates, with QBS-C0L possessing the lowest truncation error. 


46 



7— r—i— — i i-7 yn "■’ ■ ^iyj i nn,<n.ui.) 1 ' nww w ^ i 

I 


The behavior of QBS-COL is due to the different orders of approximation of the 
derivatives, i.e., O(h^) for the second derivative and O(h^) for the first 
derivative, which reduces to O(h^) for the overall approximation to the 
spatial operator. It is also interesting to note that CBS-COL gives 
identical results as the standard second order centered finite difference 
method, (CFD) , since the resulting matrices of the two methods are linearly 
related. 

It has been observed (for instance in Ref. 47) that for the solution of 
the boundary layer equations, near the outer edge (at large values of the 
normal coordinate), where the flow is nearly uniform, oscillations and/or 
overshoots in the velocity profile may occur. This nonphysical behavior, 
which can be traced to the violation of a cell Reynolds number stability 
condition could deteriorate the entire solution. Therefore an understanding 
of the behavior of the spatial approximations with respect to cell Reynolds 
number is desirable. Considered here is the constant coefficient equation 
once more with b set equal to 80 for 10, 20 and 40 intervals which correspond 
to cell Reynolds numbers of 8, 4 and 2, respectively. In Table VI these 
results are presented. Solution profiles for Rc = 10 are shown in Table VII. 
The generalized OCI schemes give uniformly monotonic solutions (they are 
constructed to do just that) , while the B-spline schemes give oscillatory 
solutions for cell Reynolds numbers approximately greater than 2. 

The behavior of the OCI schemes in the range of Reynolds number 2 to 4 
is a property of such an unrestricted (with Rc) scheme (cf., Ref. 29). As is 
true for all schemes, only in the limit as h + 0 docs one obtain the convergence 
rates predicted by the theory. Error estimates obtained by Berger et al. f 
(Ref. 29) for the generalized OCI schemes indicate there are several over- 
lapping regions that are dependent on the magnitude of Rc which varies as 

— n ^ 

h with h fixed. For p positive and large, Rc ■* 0 and one recovers an 0(h ) 

convergence rate, while for p negative and large, Rc -*• <» and one obtains an 
2 

0(h ) convergence rate. The transition, which is automatic, results in a 
second order upwind differencing formula when Rc is large. A more detailed 
discussion for the entire range of p values is presented in Ref. 29. 

By way of comparison, the exponential scheme of Il'in (Ref. 31) has been 
shown to be uniformly first order accurate (Ref. 48) while the exponential box 
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scheme due to El-Mistikawy and Werle (Ref. 35) is uniformly second order 
accurate (Ref. 36). For the problem under consideration (constant coefficient) 
both exponential schemes give the exact solution as expected. The above 
discussion points out one of the major advantages of the Q-R operator format. 
Once a program has been written within such a framework, various schemes can 
be implemented easily. Even hybrid type schemes that are evaluated pointwise 
are allowed. For example, schemes can be chosen by considerations of accuracy, 
cell Reynolds number behavior, or running time. 

For actual viscous flow problems, where moderate to high cell Reynolds 
numbers appear in uniform regions (near the outer edge of the viscous layer) , 
the generalized OCI schemes work well in practice (Ref. 49) . 

It is important to note that the cubic B-spline generalized Galerkin 
procedure is sensitive to the type of boundary condition set at x = 0. When 
collocation at x •» 0 was employed, the results were nonphysical while setting 
the first derivative to its exact value at that point gave results comparable 
to QBS-COL. The boundary condition was applied in a region of steep gradients 
so that the second order collocation approximation is not sufficient to 
prevent the erratic behavior. However, at the low cell Reynolds number range, 
the collocation boundary condition worked well as the results in Table V 
indicate. 

The second case considered is the linearized Burgers equation 


'a /ax \' 

u yv + — tanh — — u » 0 

xx i/ \ Zv ) x 


which has the exact solution 


1 '-'° nh (-ir) 


u(-co)— I , u(C0)-~ 0 


The coefficient a was set to 1/2 while v was varied; small v corresponds 
to a shock near x “ 0. Calculations were carried out in the regions 
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-5 £ x £ 0 and/or 0 £ x £ 5, with the boundary conditions set to their 
exact function values. The Burgers equation mimics a true boundary layer in 
that regions of sharp gradients (near x «* 0) corresponds to low cell Reynolds 
numbers while in the "uniform flow" region (|x| -*•=>) the cell Reynolds number 
reaches a maximum. 

Results are presented in Table VIII for v *■ 1/16 and for mesh intervals 
50, 100 and 150 for the following schemes; 0CI-G1, 0CI-G2, OCI-Swartz, CBS- 
GAL, QBS-COL and Allen's exponential scheme. The QBS-COL method again has 
the lowest truncation error for a given mesh distribution. However, the 
computation times (on a CDC 7600 machine) as presented in Table VIII indicate 
that in order to attain a given accuracy, 0CI-G2 is the most efficient scheme 
for this problem while there is little difference between OCI-G1 and CBS-GAL. 
Although QBS-COL does not fare as well, its performance car. be improved by 
storing the values of the basis functions, B^ , at the nodal pointp and not 
computing them as needed. In Table IX 0CI-G1, 0CI-G2, CBS-GAL, QBS-COL, the 
Allen exponential and the El-Mistikawy-Werle exponential schemes are compared 
for the case v = 1/24 as the cell Reynolds number is increased beyond 2. 
Again, the OCI schemes and the exponential schemes have monotonic behavior 
while the B-Spline techniques lead to oscillatory solutions. 

The results indicate that with respect to accuracy there is little to 
choose from among the generalized OCI schemes and the higher order B-spline 
procedures. However, the B-spline procedures lead to larger banded matrices 
and have cell Reynolds number stability restrictions. Note that it is 
certainly not ruled out here that at some point in the future the B-splines 
could be modified to eliminate the oscillatory behavior at high cell Reynolds 
number. Although the B-spline schemes do have the advantage of being able to 
evaluate derivatives and treat derivative boundary conditions more easily than 
the OCI methods, at present they do not seem to be as versatile as OCI at 
least on equation systems that OCI can be applied to. It is therefore 
recommended that for the problem under consideration here for the spatial 
scheme, the family of Q-R operator schemes, in particular the generalized 
OCI schemes be adopted. 

As a model problem for investigating the chosen spatial and temporal 
schemes, a coupled system of one-dimensional, nonlinear, parabolic equations 
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was constructed which exhibit typical boundary layer-like behavior, and have the 
form 

U t " u yy + A 2 u y + A 3 U + A 4 UW + A 9 W + A 6 (54) 

W t “ w yy + B z w y + B s w + B 4 UW + B a u + B e ( 55 > 


A 2 “ w , A 3 »-a 2 , A 4 «a, A 5 - -a(e + sinwt), 

A. » u»coswt(l-e ay ) + a 2 (c + sinujl) 

D 


0 2 * u > B 3 - -a , B 4 - a, 

B s * -aSsinait - ye ay (2y + /3sintut), 

B 6 ■ cucoswt ■[ yfiye ay + 8(l-e ay )j + a 2 8sinojf 
- 2 ye’ ay | ! - a(2y +/3sinajt)} 


with boundary conditions 


u(0,t) « w(0,t) « 0 


u(co ,t) » c + sinwt c > I 


w(co,t) = Ssincut 8 > 0 


The initial conditions were chosen to be the exact solution values at the 
initial time. 


The exact solutions are 


u * (c + sinwtK I - e' ay ) 


w » yy(y + /3sinwf)e ay + 8sinwt( I - e’ ay ) 
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The exact solutions possess the following properties: 

1. periodicity in time, 

2. exponential decay in space, 

3. nonlinearities and equation coupling, 

and w exhibits a form of "flow reversal". 

The Q-R operator spatial schemes were examined, including second order 
central finite differences, second order Il'in (exponential) and the fourth 
order generalized operator compact implicit method. As noted previously, all 
the above schemes require the inversion of a scalar tridiagonal matrix for a 
single equation in one unknown. 

In the computer code, Eqs. (54) and (55) were solved simultaneously, 
resulting in a 2 x 2 block tridiagonal system. Both first order fully implicit 
and second order Crank-Nicolson discretizations were included in the program. 
Appendix B presents the appropriate linearizations and describes the implemen- 
tation of the method. 

The functions u and w were chosen, as previously mentioned, to mimic the 
velocity components of a time-dependent boundary layer flow. The "u velocity 
component" is always positive, with its maximum value (which is reached at the 
outer edge) varying between e + 1 and c - 1. The "w cross flow velocity 
component", as shown in Fig. 1, has a maximum value in the interior, and 
exhibits a cross flow reversal in that w changes sign through a time cycle. 

The numerical results to be discussed are for the case shown, i.e., a » 1, 

3 “ 4.0, y " -2.0, 6 » 0.40, and e * 1.5. 

A comparison of Newton iteration with the second order noniterative time 
linearization procedure was made. The results shown in Figs. 2 and 3 verify the 
contention that the effort spent at iterating can be more effectively used in 
decreasing the step size, thus reducing both the temporal truncation and 
linearization errors. The spatial step size for the calculation shown in 
Figs. 2 and 3 was chosen (at h = .2), and the method employed was 0C1-G1, so 
that the predominant error would be due to temporal effects, e.g., linearization 
and discretization. As an auxiliary benefit, time linearization requires one 
less time level of storage, and is easier to program (cf.. Appendix B). 

Figures 4 and 5 present a comparison of the generalized 0C1 scheme with 
second order central differences. The benefits In accuracy of the higher order 
OCI scheme is evident in these figures. In fact, for the J = 40 case, the 
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temporal truncation error, which can be shown to be of order 10 - ^, dominated 
the spatial truncation error so that the OCI scheme did not attain its 
theoretical convergence rate, 
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In this report several temporal schemes have been Investigated, consis- 
tently split block Implicit, Rubin's predictor-corrector and Hopscotch. 

It has been shown that, in order to meet the requirements of three-dimensional 
unsteady boundary layer flows, the consistently split linearized block implicit 
scheme is the most versatile and efficient of the three, and was thus 
recommended. 

Two approaches to spatial approximation were described, the Q-R operator 
formulation and the B-spllnc basis function technique. Results of numerical 
experlmints indicate that both qulntic B-spline collocation and generalized 
OCI schemes performed well and that the Q-R operator formulation, in particular 
the generalized OCI schemes, are at present particularly well suited for the 
problem of time-dependent boundary layers in regard to efficiency, cell 
Reynolds number stability restrictions, and flexibility. 

The linearized block implicit temporal scheme, in conjunction with a 
generalized OCI spatial scheme, was employed to solve a system of two coupled 
nonlinear parabolic equations that exhibit "three-dimensional unsteady boundary 
layer" behavior. The results of numerical experiments indicate that the 
generalized OCI approach is viable and can be applied to viscous flow problems. 
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APPENDIX A 


Linearization Technique 

A number of techniques have been used for implicit solution of the 
following first-order nonlinear scalar equation in one dependent variable 
$(x,t) : 


<3<£/<3t =F{</>) <3G {£)/<3x 


(Al) 


Special cases of Eq. (Al) include the conservation form if F($) *» 1, ond 
quasi-lincar flow if G(<f>) “ Previous implicit methods for Eq. (Al) 
which employ nonlinear difference equations and also methods based on two- 
step predictor-corrector schemes are discussed by Ames (Ref. 50, p. 82) and 
von Rosenburg (Ref. 51), p. 56). One such method 1 b to difference nonlinear 
terms directly at the implicit time level to obtain nonlinear implicit 
difference equations; these are then solved iteratively by a procedure such 
as Newton's method. Although otherwise attractive, there may be difficulty 
with convergence in the iterative solution of the nonlinear difference 
equations, and some efficiency is sacrificed by the need for iteration. An 
implicit predictor-corrector technique has been devised by Douglas and Jones 
(Ref. 52) which is applicable to the quaslllnear case (G » $) of Eq. (Al). 

The first step of their procedure is to linearize the equation by evaluating 
the nonlinear coefficient as F($ n ) and to predict values of using either 

the backward difference or the Crank-Nicolson scheme. Values for 4> n+1 are 
then computed in a similar .manner using F($ n+ *^) and the Crank-Nicolson scheme. 
Gourlay and Morris (Ref. 53) have also proposed implicit predictor-corrector 
techniques which can be applied to Eq. (Al). In the conservative case (F ■ 1), 
their technique is to define G($) by the relation G($) » $G(<J) when such a 
definition exists, and to evaluate G($ n+ *) using values for ^ 1 computed by 
an explicit predictor scheme. With G thereby known at the implicit time level, 
the equation can be treated as linear and corrected values of are computed 

by the Crank-Nicolson scheme. 

A technique is described here for deriving linear implicit difference 
approximations for nonlinear differential equations. The technique is based 
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on an expansion of nonlinear implicit teres about the solution at the known 
time level, t n , and leads to a one-step, two-level scheme which, being linear 
in unknown (implicit) quantities, can be solved efficiently without iteration. 
This idea was applied by Richtmyer and Morton (Ref. 16, p. 203) to a scalar 
nonlinear diffusion equation. Here, the technique is developed for problems 
governed by l nonlinear equations in t dependent variables which are functions 
of time and space coordinates. The technique will be described for the three- 
dimensional, unsteady equations. 

The solution domain is discretized by grid points having equal spacings 

12 3 12 3 

in the computational coordinates. Ay , Ay and Ay in the y , y and y 

directions, respectively, and an arbitrary time step. At. The subscripts i, j, 

12 3 

k and superscript n are grid point Indices associated with y , y , y and t, 

respectively, and thus ^ ^ denotes $(y^, yj, y^» t n ). It is assumed that 

the solution is known at the n level, t n , and is desired at the (n+1) level, 

t n+ *. At the risk of an occasional ambiguity, one or more of the subscripts 

is frequently omitted, so that $ n is equivalent to , 

*»J 

The numerical method employed is quite general and is formally derived for 
systems of governing equations which have the following fora: 


= (<£) + S (£) 


(A2) 


where $ is a column vector containing l dependent variables, H and S are 
column vector functions of if, and .2) is a column vector whose elements are 
spatial differential operators which may be multidimensional. The generality 
of Eq. (A2) allows the method to be developed concisely and permits various 
extensions and modifications (o.g., noncartcsian coordinate systems, turbulence 
models) to be made more or less routinely. It should be emphasized, however, 
that the Jacobian 311/3$ must usually be nonsingular if the ADI techniques as 
applied to Eq. (A2) are to be valid. A necessary condition is that each 
dependent variable appear in one or more of the governing equations as a time 
derivative. An exception would occur if for Instance, a variable having no 
time derivative also appeared in only one equation, so that this equation could 
be decoupled from the remaining equations and solved a posteriori by an alter- 
nate method. 
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The linearized difference approximation is derived from the following 
implicit time-difference replacement of Eq. (A2): 

(H n + '-H n )/AI=/3[2»(^ n+, )+S n + , ]+(l-/3)[2>( <#> n )^S n ] (A3) 

where, for example, H n+1 i . The form of 2> and the spatial differ- 

encing are as yet unspecified. A parameter 6(0 5 6-1) has been introduced 
so as to permit a variable centering of the scheme in time. Equation (A3) 
produces a backward difference formulation for 6 “ 1 and a Crank-Nicolson 
formulation for 6 “ 1/2. 

The linearization is performed by a two-step process of expansion about 
the known time level t n and subsequent approximation of the quantity 
(34>/3t) n At, which arises from chain rule differentiation, by ($ n+ ^ - $ n ). The 
result is 

H n+I = H n +(dH/d£) n (<£ n + '-<£ n ) + 0 ( Af) Z (A4a) 

s n+ ' = s n +(as/^) n (* n+, -<n+0(At) 2 (A4b) 

2>(* n+ ') = %(<}> n )+(dZ/d4>){<t> n + '-<t> n )+0{M) Z (A4c) 

The matrices 3H/3if> and 3S/3$ are standard Jacobians whose elements are defined, 

for example, by (3H/3$)^ r = 311^/3$^. The operator elements of the matrix 

3 ^/ 31 *, are similarly ordered, i.e., (3.2>/3$) Z 3^ / 3<^ ; however, the 

cjr q r 

intended meaning of the operator elements requires some clarification* For 
the q 11 * 1 row, the operation (3.2>^/3(|>) n ($ n+1 - $ n ) is understood to mean that 
(3/3t2) [#(x,y,z,t)]) n At is computed and that all occurrences of (3$ r /3t) n 
arising from chain rule differentiation arc replaced by (<}>” +1 - $")/At. 

After linearization as in Eqs. (A4) , Eq. (A3) becomes the following linear 
implicit time-differenced scheme: 
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{aH r 7<#W n+, -£ n )/Af =J2>{«^ n > +s ft + p (d%/d$ + ds n /d<}‘U<t> n+ i-<i> n ) (as) 

Although H° + * is linearized to second order In Eq. (A4) , the division by At 
In Eq. (A3) Introduces an error term of order At. A technique for maintaining 
formal second-order accuracy in the presence of nonlinear time derivatives is 
discussed by McDonald and Briley (Ref. 7), however, a three-level scheme 
results. Second-order temporal accuracy can also be obtained (for B ■ 1/2) by 
a change in dependent variable to $ = H(d>) , provided this is convenient, since 
the nonlinear time derivative is then eliminated. The temporal accuracy 
is independent of the spatial accuracy. 

On examination, it can be seen that Eq. (A5) is linear in the quantity 
($ n+ ^ - $ n ) and that all other quantities are either known or evaluated at 
the n level. Computationally, it is convenient to solve Eq. (A5) for 
($ n+ * - $ n ) rather than $ n+ ^. This both simplifies Eq. (A5) and reduces 
roundoff errors, since it is presumably better to compute a small 0(At) change 
in an 0(1) quantity than the quantity itself. To simplify the notation, a 
new dependent variable defined by 

— <£ n (A6) 

is introduced, and thus ifi n+7 “ - $ n , and i^ n ■ 0. It is also convenient 

to rewrite Eq. (A5) in the following simplified form: 

(A + At^)'/' n + l = At [- 2 > (£ n )+S n ] (A7a) 

where the following symbols have been introduced to simplify the notation: 

As dH n /d<f> -/?At(c?S n A)<£) (A7b) 

1~ -p{d (A7c) 

It is noted that J!(\ p) is a linear transformation and thus J{0) “ 0. Further- 
more if $) is linear, then m -B^>(i^). 



Spatial differencing of Eq. (A7a) is accomplished simply by replacing 

i 2 i i 

derivative operators such as 3/3y , 3 / 3y 3y by corresponding finite 
difference operators, D^, nenceforth, it is assumed that J2) and 1 have 

been discretized in this manner, unless otherwise noted. 

Before proceeding, some general observations seem appropriate. The 
foregoing linearization technique assumes only Taylor expandability, an assump- 
tion already implicit in the use of a finite difference method. The governing 
equations and boundary conditions are addressed directly as a system of coupled 
nonlinear equations which collectively determine the solution. The approach 
thus seems more natural than that of making ad hoc linearization and decoupling 
approximations, as is often done in applying implicit schemes to coupled 
and/or nonlinear partial differential equations. With the present approach, 
it is not necessary to associate each governing equation and boundary condition 
with a particular dependent variable and then to identify various "nonlinear 
coefficients" and "coupling terms" which must then be treated by lagging, 
predictor-corrector techniques, or iteration. The Taylor expansion procedure 
is analogous to that used in the generalized Newton-Raphson or quasi- 
linearization methods for iterative solution of nonlinear systems by expansion 
about a known current guess at the solution (e.g.. Bellman & Kalaba, Ref. 54). 
However, the concept of expanding about the previous time level apparently 
had not been employed to produce a noniterative implicit time-dependent scheme 
for coupled equations, wherein nonlinear terms are approximated to a level of 
accuracy commensurate with that of the time differencing. The linearization 
technique also permits the implicit treatment of coupled nonlinear boundary 
conditions, such as stagnation pressure and enthalpy at subsonic inlet 
boundaries, and in practice, this latter feature was found to be crucial to 
the stability of the overall method (Ref. 17). 
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APPENDIX B 

Discretization of Model Equations 

Herein we shall consider the system of model equations discussed in the 
section on numerical experiments namely, 


u yy + A 2 u y + A s u + A 4 UW + A s w + A 6 


(Bl) 


where 


v, ■ w ,. A 3 ■ -a , A 4 «a, A 9 - -a(e + sinwt), 

A ■ o)coswt{l-e" ay ) + a 2 U + sinwt) , 

v 


and 


w t “ w yy + B 2 w y + B a w + B 4 UW + B s u + B c 


(B2) 


where 


B 4 - a, 


B 2 ■ U, B, - -a , 

B 8 ■ -aSsinwt -ye ay (2y +/3sinwt), 

B g • wcoswl | yfiye ay + 8(l-e ay )| + a 2 8sinu>t 
- Sye^^l-alSy+^sInwt)} 

Since nonlinear terms wu^, uw^ and uw appear in Eqs. (Bl) and (B2) , 
linearization is required. Two types of linearization procedures are described, 
time linearization and Newton iteration (quasilinearization). The application 
of these methods within the Q-R operator framework is also demonstrated. 


59 


4 - 



nrur t xi 177J3LE»\: 


r' y '-' v 

&\- 

J.S 


. T . . .■ , ; r ■ 




Time Linearization 

The nonlinear terms uv and wu are linearized by the method described in 
Appendix A to yield 


(w V n+/3 - w n (/3u y n+ ' + d-/9)u" ) + / 3 ujw" 4 ' -/ 3 w n Uy n 


(uw) n+/? - /3wV + ' + /3u n w nt| + (i-2/3)w n u n 

where 6 - 1/2 corresponds to a Crank-Nicolson scheme and 0 - 1 to a fully 
implicit method. 

Substituting Eqs, (B3) and (B4) into (Bl) we obtain 

(u n4 '-u n )/A» . -^-[/3u n+l + (l-/3)u n ] 

+ A 5 [/3w n u n *' + /W 1+ ' + (|-2/?)u n w n ] 

+ A,[ ^w n -' + (|-£) w n ] + a 7 + /3u n w n*i . ^ n 


(B3) 


(B4) 


which, with some rearrangement, reduces to 

I 


M X 0 - Ay 2 ( A s w n )] -R}u n+I - Ay 2 Q{u y n + A 3 u n + A 6 } w 

I . z ( (1-2/9) \ 

~P + ** ( A , + A 5 w j 


(1-/9) 

+ — Q - R 


0-/3) 


+ Ay z Q( — - — a w" + - u, „ 

7 P 6 W p -W U y 


where 


X « At /Ay 2 

and Q and R are evaluated at the n th time level. 


(B5) 
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Similarly, Eq. (B2) reduces to 




{°[t« - B 5 u")] - r}w"*' - 4y*0 


’[ T? + B 5 * ^■V") 


+ b 5 » + B, 


n+l 


( I ~P) 

+ __ R 


w 


+ 4y '° ^ V" + - u'S" 


(B6) 


Equations (B5) and (B6) are solved as a coupled (2 x 2 block) system for u n+1 


j n+l 
and w 


Newton Iteration 

Alternatively one can linearize about a previous iteration instead of the 
previous time level, to yield for the nonlinear term (wu ) n+6 for example, 


(w *V * [ + /3<i-/3)w"]u y n4 ' + [/3 2 Uy * + /3(i-/3) Uy n ] w n 

+ [(l-/3) 2 w n u y n -/3Vu y *] 


(B7) 


/ 




/ 


Substituting Eq. (B7) in Eq. (Bl) the quasilinear approximation becomes 

(u n4 ' -u n )/At * [/3u n41 + d-/3)u n ] 

+ A 5 [(/3w* + (i-/3) w n )/3 u n41 +(£u*+(|-/3)u n )/3w rHI 
+ ( ! -/3) 2 u n w n - /? 2 U V* + A 6 [/3 w" 4 ' + (l-/3)w n l 

i ^ 

+ A 7 + [£ Uy * +(l-/3)u y n ]/3w n41 -[/3u* + n-/3)u y n ]/3w* 
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which after some manipulation reduces to 


n+l 


{o[~- Ay z (A s (/3w* + (l-/3)w n ))] -r} u 
- Ay 2 o{A 3 (/3u* + (l-/3)u n ) +A 6 + (/3u y # +(l-/3)u y n )}w 


{o [_l + ^^ (V ,.^ )] + ^ r } u 


(I-/?) 


(1-/3) 


0 


+ Ay 2 o{-^- + ^-^-A 6 w n - /3A 5 u*w*- (/3u y # +(|-y9)u y n )w*} 


(B8) 


Here Q and R are assumed to be evaluated at (n+g) . A similar expression 
can be obtained for Eq. (B2) , viz., 


{o[ — - Ay 2 (B 5 (/3u^-(l-/3)u n ))]-R}w n+l 
-Ay z Q {b 5 (£w* +( l-£)w n ) +B 6 +(/3vv y * +d-/3)w y )}u n+l 

'Hip + Ay'^tBgd-^Ju")] + r} w n 

+ Ay z o{^- + B 6 u"-/3B 5 u tt w*- (£w* + (l-/3)w y )u*} 


(B9) 


Equations (B8) and (B9) are solved as a coupled system for u and w , the 
latest lterants . The coefficients and the Q and R operators are updated, and 
the system of equations is solved again. The process is continued until con- 
vergence is attained. The process involves : 

1. updating the coefficients Ay - Ay and By ~ By, 

2. updating the 0 and R operators, 

3. solving the resulting 2x2 block system, 

4. repeating steps 1-3 until convergence is attained. 


t 
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Note that Eqs. (B5), (B6), (B8) and (B9) can easily be transformed to increment 
fora. However, the results in this report were obtained for the form of the 
equations given above. 
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TABLE II. - OPERATOR COEFFICIENTS FOR GENERALIZED 

OPERATOR COMPACT IMPLICIT SCHEME OCI-G2 

qj“ - e + [ p, - 3] RCj +[p 2 ]r c j 2 
Qj C " 60 +[|0p,]r Cj + [ P 3 ] Rc f + [ 7 "j«|P'»] Rc j* 
qj + -6 + [ P , + 3 ]rc j +[p t + p 2 ]Rc J 2 +[pjRc J 3 

where 

p, • 3 , p 2 - 0 , p 3 » mox [ 7Tj | TTg ] 

' 7r i " (t J+i + V' )p 2 + T J + ' P « + W \ w * " 5 - 2p 2 + ( t r 2 - | )p 2 -3(T, t| + o- 2 ) + w z 

0 o-|>0 _ 0 2p,-<r 2 >0 

-g- <r, z {IO “ Tj +I * Tj_, ) (Tj <0 2 (2p,-cr 2 ) 2 /8 2p,-cr 2 <0 

"i S'. •3r | „-T J ., + ,O t2 h T) .,(iL.) 

p * - Tt l + r j*i] ’''s * ’’'i + + 2T i-i( 2 + h T~') p t 

with h sufficiently snail so that 

IObj-bj.,-bj +l > 0 and hCj 41 /bj 4 ,<2 for j»2,***,J and CjSO 

•here 

r J-i “ b j-i /b J * T J*i * b J*i /b J 0nd Rc J " hb J 
fj, rj, r| given in TABLE I 
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TABLE III. - OPERATOR COEFFICIENTS FOR GENERALIZED 
OPERATOR COMPACT IMPLICIT SCHEME OCI-G1 


qf ,c,+ ■ °o" ,c,+ + Q r ,c,+ + ° 2 " ,c, ‘ f + Q i ,c '* 

Qq -60 » Qq ■ Qq * ® 


,C 


'I 10 


c /r_ p_ p \ 

" Q o\ 12 150 75 / 


. pC P ~ 2P ' \ 

Q t 10 °°\ 12 150 75 / 


O’ « min (/ l,,/x 2 ) 


Qj - a< 2 - ZQ- - -g-(- 2oJ> C + 0 *p* - II Q ~p~) 

Q z " Q z + T (0 .V - ZQ iP~ - 2Q ,V> 


°3 ’ l w 3 + [ 3Q 2 >" “ 0 2 V + ° 2 V]}/ (I + P* / P C) 


°3 ’ 


\ 



T/.. 


TABLE III. - CONTINUED 


where 


o» 2 ■ min(/x 3 , /i 4 ) 


w, • 0 


- \2 


M, -Co-1 


H-z " -2Q,V) 


5 r 6 _c n c/ Z 5 ' 13 _ 3 ♦yi 5 

^3 ’ 240q l 5 °l ’ °0 V 2 IOO ^ 1 00 ^ /i 


I 

^4 " 240 qP C 


(q y - q ,>’ -°iV)/[/» e ’ iV ( ^ ,4+/5 ’ ) ] 


+ 4(/> 4 + + Zp C )Q, - 0 


+ y[(- 2 Q ^ C + 0 |V " 


+ -^c (0^ C “ 20," /)- - 20, V)] 


with h sufficiently small so that 


l0b ] • b H • b j*. > 0 


r j* r j’ r J coc f f ielents are given in TABLE I 
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TABLE IV. - OPERATOR COEFFICIENTS FOR EL-M1STIKAWY WERLE 
EXPONENTIAL BOX SCHEME 



P~ ' -T ( P\-\ + P\ ] » P* * + 

P\ ' hbj 


72 




































■■■■•" a ■ ' Q 


TABLE VI. - COMPARISON OF SPATIAL ERRORS FOR 
CONSTANT COEFFICIENT PROBLEM 


Scheme 

Jmax 

Rc 

Max error 

Comment 

CBS-CAL 1 

40 

2 

-.158 E+01 

Oscillatory 

QBS-COL 

40 

2 

-.484 E-03 

Oscillatory 

0CI-G1 

40 

2 

-.162 E-01 

Mortotone 

OCI-G2 

40 

2 

-.863 E-02 

Monotone 

CBS-CAL 1 

20 

4 

-.145 E+01 

Oscillatory 

QBS-COL 

20 

4 

.612 E-01 

Oscillatory 

OCI-G1 

20 

4 

-.586 E-01 

Monotone 

OCI-G2 

20 

4 

-.127 E-01 

Monotone 

CBS-GAL 

10 

8 

-.136 E+01 

Oscillatory 

QBS-COL 

10 

8 

.288 E+02 

Oscillatory 

OCI-Gl 

10 

8 

-.561 E-01 

Monotone 

OCI-G2 

10 

8 

-.425 E-02 

Monotone 


^Collocation at x=0 used. 


• \ 
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TABLE VII. - SOLUTION PROFILES FOR CONSTANT 
COEFFICIENT PROBLEM - Rc - 10 


X 

Exact Solution 

Computed Solution 



(QBS-C0L) 

(CBS-GAL) 1 

(OCI-G2) 

0.0 

0 . 

0 . 

0 . 

0 . 

0.1 

.999955 E+00 

-.374802 E+00 

.118628 E+01 

.997532 E+00 

0.2 

.999999 E+00 

.304662 E+01 

.965301 E+00 

.999994 E+00 

0.3 

.100000 E+01 

.311375 E-01 

.100646 E+01 

.999999 E+00 

0.4 

.100000 E+01 

.251958 E+01 

.998796 E+00 

.100000 E+01 

0.5 

.100000 E+01 

.476597 E+00 

.100022 E+01 

.100000 E+01 

0.6 

.100030 E+01 

.215299 E+01 

.999958 E+00 

.100000 E+01 

0.7 

.100000 E+01 

.778851 E+00 

.100001 E+01 

.100000 E+01 

0.8 

.100000 E+01 

.189431 E+01 

.999999 E+00 

.100000 E+01 

0.9 

.100000 E+01 

.107440 E+01 

.100000 E+01 

.100000 E+01 

1.0 

.100000 E+01 

.100000 E+01 

.100000 E+01 

.100000 E+01 


^Exact derivative at x ■ 0 set. 


75 


























f 




t- 
i : 


i m i .i 

»T"waw i ." ^ , "'j | ."» « 1 »■■.■- ■■ 

' ' ■ • a ■ 


TABLE IX. - COMPARISON OF SPATIAL ERRORS FOR 
BURGERS EQUATION - V - 1/24 


Scheme 

Jmax 

Rc max 

Max error 

Comment 

Allen 

25 

2.4 

.182 E-01 

Monotone 

El-Mistikawy 

Werle 

25 

2.4 

-.213 E-01 

Monotone 

CBS-GAL 

25 

2.4 

.242 E-02 

Oscillations 

QBS-COL 

25 

2.4 

.243 E-02 

Oscillations 

OCI-G1 

25 

2.4 

.438 E-02 

Monotone 

OCI-G2 

25 

2.4 

.372 E-03 

Monotone 

Allen 

50 

1.2 

-.814 E-02 

Monotone 

El-Mistikawy 

Werle 

50 

1.2 

-.758 E-02 

Monotone 

CBS-GAL 

50 

1.2 

-.800 E-03 

Monotone 

QBS-COL 

50 

1.2 

.392 E-03 

Monotone 

OCI-G1 

50 

1.2 

.123 E-02 

Monotone 

OCI-G2 

50 

1.2 

-.279 E-03 

Monotone 
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Figure 1. - Cross flow velocity profiles as a function of time. 
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Figure 2. - Maximum 
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